#!/usr/bin/perl
# 2012/11/30 Estimate HiD resolution

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use File::Basename;
use Carp qw(croak);
$| = 0;

if(@ARGV != 2 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file]\n";
}

my %opt;
getopts("i:", \%opt);
my $FILE_map = $opt{i};
my ($map_name, $map_dir, $map_ext) = &fileparse($FILE_map, '\..*');

my $PROGRAM_DIR = '/cscb/noma_nfs/Program/HiD/';
my $FILE_HindIII_sites = '/cscb/noma_nfs/Project/2012-11-15_HiD_1%4%/out/HindIII_sites.txt';
my $PROGRAM_conversionTalbe = $PROGRAM_DIR . 'Make_conversionTable_HindIII2resolution.pl';
my $FILE_unAlignable_hinID = '/cscb/noma_nfs/Project/2012-11-15_HiD_1%4%/out/2012-12-06_un_alignable_fragment_hinID_list.txt';


my @CHECK_RESOLUTIONS = qw(2000 3000 4000 5000 6000 7000 8000 9000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 70000 80000 90000);
my $TMP_conversionTable = 'tmp_' . $map_name . '_conversion_table.txt';



#---------------------------------------
# read un-alignable hinID list
#---------------------------------------
my %unAlign;
my $fh_un = IO::File->new($FILE_unAlignable_hinID) or die "cannot open $FILE_unAlignable_hinID: $!";
while($_ = $fh_un->getline()){
	s/\r?\n//;
	$unAlign{$_} = 1;
}
$fh_un->close();



print "resolution\tcombination\tcis num\t%\ttrans num\t%\n";
foreach my $resolution (@CHECK_RESOLUTIONS){
	my ($total, $cis, $trans) = &Calc($resolution);
	printf "$resolution\t$total\t$cis\t%.2f\t$trans\t%.2f\n", $cis/$total*100, $trans/$total*100;
}

unlink($TMP_conversionTable);

sub Calc{
	my ($resolution) = @_;

	#---------------------------------------
	# make conversion table
	#---------------------------------------
	system("perl $PROGRAM_conversionTalbe -i $FILE_HindIII_sites -o $TMP_conversionTable -b $resolution") == 0 or die "conversion failed: $!";

	#---------------------------------------
	# read conversion table
	#---------------------------------------
	my %table;
	my @binLists;
	my $fh_conv = IO::File->new($TMP_conversionTable) or die "cannot open $TMP_conversionTable: $!";
	while($_ = $fh_conv->getline()){
		s/\r?\n//;
		my ($id, @hin) = split /\t/;

		my $total_ok_hin = 0;
		foreach my $h(@hin){
			unless(exists $unAlign{$h}){
				$table{$h} = $id;
				$total_ok_hin++;
			}
		}
		if($total_ok_hin > 0){
			push @binLists, $id;
		}
	}
	$fh_conv->close();

	#---------------------------------------
	# correct information
	#---------------------------------------
	my %data;
	my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
	while($_ = $fh_map->getline()){
		s/\r?\n//;
		my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

		unless(exists $table{$hinID1}){
			next;
		}
		unless(exists $table{$hinID2}){
			next;
		}

		my $bin1 = $table{$hinID1};
		my $bin2 = $table{$hinID2};

		my $hinNum1 = substr($hinID1, 0, -1);
		my $hinNum2 = substr($hinID2, 0, -1);

		if($hinNum1 < $hinNum2){
			$data{$bin1}{$bin2}++;
		}else{
			$data{$bin2}{$bin1}++;
		}
	}
	$fh_map->close();


	#---------------------------------------
	# output
	#---------------------------------------
	my $TOTAL = 0;
	my $CIS = 0;
	my $TRANS = 0;
	for(my $i = 0; $i < @binLists; $i++){
		my $bin1 = $binLists[$i];
		my $chr1;
		if($bin1 =~ m/bin_\d+\|(chr\d)/){
			$chr1 = $1;
		}
		for(my $j = 0; $j < @binLists; $j++){
			my $bin2 = $binLists[$j];
			my $chr2;
			if($bin2 =~ m/bin_\d+\|(chr\d)/){
				$chr2 = $1;
			}
			my $value = 0;
			if($i < $j){
				$value = exists $data{$bin1}{$bin2} ? $data{$bin1}{$bin2} : 0;
			}else{
				$value = exists $data{$bin2}{$bin1} ? $data{$bin2}{$bin1} : 0;
			}

			$TOTAL++;
			if($value != 0){
				if($chr1 eq $chr2){
					$CIS++;
				}else{
					$TRANS++;
				}
			}
		}
	}

	return ($TOTAL, $CIS, $TRANS);
}



